%% Load
clear;
close;

[name path] = uigetfile({'*.GTB'});
Head = file_Load(name, path, 'Head');
Data = file_Load(name, path, 'Data');
Data = util_DataReference(Data, Head.RefeChan);

Fig_ScrollData(Head, Data);

%% Filter
fs = Head.SampRate;
Wn = [10 16];
win = 'hann';
mul = 5;

for c = 1 : Head.ChanNum
    [Data(c,:) b] = FirFilter(Data(c,:),fs,Wn,win,mul);
end
figure;freqz(b);

Fig_ScrollData(Head, Data);





%% 1st criterion: Thresholding P < 0.01
fs     = Head.SampRate;
chan   = util_GetIndex(Head.ChanLabel, {'C3'});
X      = Data(chan,:);
Z      = zscore(abs(hilbert(X)));
N2     = util_GetIndex(Head.Stage.Label, {'N2'});
epch   = Head.Stage.Series(:,N2) == N2;
stime  = Head.StartDate + Head.StartTime;

peaks = GetPeak(Z,1);
candi = peaks(Z(peaks) > 2.328); % One-tailed P = 0.05 -> 1.645, P = 0.01 -> 2.328

spindle = [];
for e = 1 : length(epch)
    if epch(e)
        etime = Head.Stage.Time(e);
        sidx = util_GetOffIndex(stime,fs,etime);
        eidx = sidx+6000-1;
        [~,idx] = find(bitand(candi >= sidx, candi <= eidx));
        spindle = [spindle candi(idx)];
    end
end

%% 2nd criterion: Spindle Duration
z   = Z(spindle);
cut   = 0.5;          % spindle duration -> the half of the peak
dur   = 300;          % dur > 400 ms

rej = [];
for s = 1 : length(spindle)
    l = 0;
    while Z(spindle(s)-l) >= (z(s)*cut)
        l = l + 1;
    end
    r = 0;
    while Z(spindle(s)+r) >= (z(s)*cut)
        r = r + 1;
    end

    sample = Z(spindle(s)-l:spindle(s)+r);
    n_sample = length(sample);
    
    % 2nd criterion
    if (n_sample*(1000/fs)) < dur
        rej = [rej s];
    end
    
end
spindle(rej) = [];

%% 3rd criterion: Shape
z   = Z(spindle);
alpha = 0.01;         % P val < 0.05

rej = [];
for s = 1 : length(spindle)
    l = 0;
    while Z(spindle(s)-l) >= (z(s)*cut)
        l = l + 1;
    end
    r = 0;
    while Z(spindle(s)+r) >= (z(s)*cut)
        r = r + 1;
    end

    sample = Z(spindle(s)-l:spindle(s)+r);
    n_sample = length(sample);
    
    % 3rd criterion
    win = triang(n_sample);
    [R P] = corrcoef(sample,win);
    if P(2) > alpha
        rej = [rej s];
    end
    
end
spindle(rej) = [];


%% remove overlap
s = 1;
while s < length(spindle)
    if (spindle(s+1) - spindle(s)) < (2*dur/(1000/fs))
        spindle(s+1) = [];
    else
        s = s + 1;
    end
end


%% Export to excel file

t = util_GetOffTime(Head.StartDate+Head.StartTime,Head.SampRate,spindle);
out = datestr(t,'HH:MM:SS.FFF');
% xlswrite('spindle.xlsx',out);


%% 
z   = Z(spindle);


l = 0;
while Z(spindle(s)-l) >= (z(s)*cut)
    l = l + 1;
end
r = 0;
while Z(spindle(s)+r) >= (z(s)*cut)
    r = r + 1;
end

sample = Z(spindle(s)-l:spindle(s)+r);
figure;plot(sample);
figure;plot(X(spindle(s)-l:spindle(s)+r));set(gca,'ydir','reverse');



%%
stime          = Head.StartDate + Head.StartTime;
Event.Label    = 'Spindle';
Event.Time     = util_GetOffTime(stime,Head.SampRate,spindle);
Event.Channel  = 'C3';

Head.Event = Event;
Fig_ScrollData(Head, Data);